FUB-HEP 19/92 
HLRZ Preprint 77/92 



Critical Exponents of the Classical 
3D Heisenberg Model: 

A Single-Cluster Monte Carlo Studyfl 

Christian Holm 1 and Wolfhard Janke 2 

1 Institut fur Theoretische Physik, Freie Universitat Berlin 
Arnimallee 14, D-1000 Berlin 33, Germany 

2 

Hochstleistungsrechenzentrum, Forschungszentrum Jiilich 
Postfach 1913, D-5170 Jiilich, Germany 

Abstract 

We have simulated the three-dimensional Heisenberg model on sim- 
ple cubic lattices, using the single-cluster Monte Carlo update algo- 
rithm. The expected pronounced reduction of critical slowing down 
at the phase transition is verified. This allows simulations on sig- 
nificantly larger lattices than in previous studies and consequently a 
better control over systematic errors. In one set of simulations we 
employ the usual finite-size scaling methods to compute the critical 
exponents v, a, /3, 7, 77 from a few measurements in the vicinity of the 
critical point, making extensive use of histogram reweighting and opti- 
mization techniques. In another set of simulations we report measure- 
ments of improved estimators for the spatial correlation length and the 
susceptibility in the high-temperature phase, obtained on lattices with 
up to 100 3 spins. This enables us to compute independent estimates 
of v and 7 from power-law fits of their critical divergencies. 



*Work supported in part by Deutsche Forschungsgemeinschaft under grant K1256. 



1 Introduction 



Critical exponents are the distinguishing parameters characterizing continu- 
ous phase transitions. Most theoretical estimates of their values have been 
calculated along three different routes. First, assuming universality, from (re- 
summed) perturbation expansions of generic field theoretical models. Second, 
from high-temperature series expansions of lattice models, and third, from 
Monte Carlo (MC) simulations of the associated Boltzmann distributions. 
All these approaches are based on resummation or extrapolation techniques 
whose systematic errors are difficult to control. To gain confidence in the 
numerical values of the exponents it is therefore quite important to have 
several independent approaches for cross-checks available. 

Until a few years ago the precision of Monte Carlo estimates was com- 
paratively poor, since this approach was plagued not only by systematic but 
also by statistical errors. The problem is that MC algorithms based on lo- 
cal update procedures are severely hampered near criticality by extremely 
long autocorrelation times (so-called critical slowing down) which reduce the 
effective statistics and thus increase the statistical errors considerably [[!]]. 
The recent development of global update algorithms that overcome the 
problem of critical slowing down is a major step forward. It is this algo- 
rithmic improvement combined with the higher speed of modern computers 
that makes systematic tests of finite-size scaling (FSS) predictions |§ and 
a reliable computation of critical exponents much more feasible than a few 
years ago. 

For spin systems cluster algorithms j|, [| turned out to be particularly 
successful. Recent applications of the multiple and single || cluster vari- 
ants to two-dimensional (2D) Ising [|, [?|, ||, XY || [TI], [II], [12|], Heisenberg 



[|13| , and other 0(n) [[nj models have demonstrated that both variants 
are equally good tools to simulate these 2D systems. In three dimensions 
(3D), however, extensive tests for the Ising || |7], and XY [IS] models 



clearly showed that the single-cluster variant is the superior algorithm. In 
our study of the 3D Heisenberg model we have therefore chosen the single- 



cluster update algorithm |19j. The set-up of our simulations is described in 
some detail in Sec. 2. 

There are of course many sources for a comparison with simulations based 



on the standard local heatbath p0[ or Metropolis [El], [22|] algorithm. In a 



recent series of papers Peczak et al. p3| |24|] used the latter algorithm com 
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bined with histogram reweighting techniques |]2"5| , |26| , |2"7 , [28| 1 to estimate the 
critical coupling and the critical exponents of the 3D Heisenberg model from 
a FSS analysis of simulations on simple cubic lattices [2~3"| and investigated 
also its (exponential) autocorrelation time r |[24|| . Their result, r oc L 2 , 
with dynamical critical exponent z = 1.94(6), shows that r is rapidly in- 
creasing with the linear lattice size L and explains why they could not go to 
systems larger than 24 3 . Since, as expected, the autocorrelation time for the 
single-cluster algorithm turns out to be almost independent of L we could 
study much larger lattices of size up to 48 3 in reasonable computer time. Our 
results described in Sec. 3 provide evidence that the asymptotic FSS region is 
indeed reached quite early or, in other words, that the amplitudes of correc- 
tion terms are very small. The present study thus significantly reduces the 
danger of systematic errors in the MC estimates of critical exponents from 
FSS analyses. 

In Sec. 4 we report another set of high-precision simulations, done this 
time in the high-temperature phase. There we use variance reduced "cluster 



estimators" |L|] for the spatial correlation length, £, and the susceptibility, x- 
By going to very large systems with up to 100 3 spins we have made sure that 
these data have only negligible finite-size corrections well below the statistical 
errors of about 0.2%. Least-square fits of the critical divergencies of £ and 
X to the well-known power-laws yield independent estimates for the critical 
coupling and the critical exponents v and 7. Obviously this is another useful 
check on residual systematic errors. Finally, in Sec. 5, we briefly discuss and 
summarize our main results. 



2 The model and simulation techniques 

Let us start with a brief description of the model and some remarks on the 
simulation techniques we have used. The partition function of the Heisenberg 
model is given by 

(1) 

where j3 = J/ksT is the (reduced) inverse temperature and the energy is 

E = J2[l-s(x)s(x + i)]. (2) 

X,i 



x 



d(f)(x)d cos 6(x) 

47T 
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Here s = (sin # cos 0, sin ^ sin 0, cos ^) are three-dimensional unit spins at the 
sites a; of a simple cubic lattice of size V = L 3 , and i are unit steps in the 
three coordinate directions. We always employ periodic boundary conditions. 



2.1 Algorithm 

For the closely related 3D Ising and XY models it has been shown || [16|, [17 



[18] that the single-cluster update is the fastest MC algorithm available. We 
have therefore chosen this variant for our study of the 3D Heisenberg model. 
One update in the single-cluster variant consists of choosing a random mirror 
plane and a random site, which is the starting point for growing a cluster of 
reflected spins. The size and shape of the cluster is controlled by a Metropolis 
like accept/reject criterion satisfying detailed balance [|, |J. Compared with 
the multiple-cluster algorithm this variant is technically somewhat simpler 
to implement and, more importantly, in three dimensions numerically more 
efficient. The reason is that, on the average, larger clusters are moved. 

To test the performance of the algorithm we have recorded autocorrelation 
functions, A(k) = p(k)/p(0), with 

p(k) = (OA +fe > - (O,) 2 , (3) 

and Oi denoting the i-th measurement of an observable. We have focussed 
on the integrated autocorrelation time, r = | + Y^T=i A(k), which describes 

the enhancement of the statistical error, e = ^Ja 2 /N V2t, for the mean value 
over a sample of N correlated measurements of variance a 2 . The infinite sum 



was always self-consistent ly cut off at /c max ~ 6r |29], |30], |3T . 

Recall that for the single-cluster update some care is necessary in defining 
the unit of time, since in each update step only a relatively small fraction 
of the spins is moved, depending on temperature and lattice size. More 
precisely, our results show that near j3 c and for all lattice sizes, the average 
cluster size, (\C\), is proportional to the susceptibility, 

(|C|>«0.75x (4) 

(with roughly the same constant as in two dimensions [14[]), where x — x/ft — 
V(m 2 ), m = yJ2xs( x )- At p c , the susceptibility behaves like x °t L 1 ^ = 
L 2 ^ 11 (with very small rj pa 0.04), so that with increasing lattice size the 
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fraction of moved spins in each update step decreases like {\C\)/V oc L~( 1+v \ 
i.e., roughly oc 1/L. Since the CPU time needed for a single-cluster update 
is approximately proportional to the number of moved spins, it is convenient 
to use N = V I (|C|) oc L 1+r > single update steps as unit of time. This is then 
directly comparable with other schemes that attempt moves for all spins in 
one update step. All our autocorrelation times refer to this unit of time 
(Metropolis-equivalent unit) which can always be achieved by a rescaling of 
the time variable. 

For the susceptibility we typically find r « 1.5 — 2.0. The values of r for 
each simulation are given in Table 1. Already our rough estimate of r shows 
that for large system sizes the single-cluster update clearly outperforms the 
Metropolis algorithm, for which the exponential autocorrelation time To has 
recently been determined |24| to be r = aL z , with amplitude^ a « 3.76 and 
dynamical critical exponent z = 1.94(6). For our largest lattice size L=48 
this implies a reduction of the autocorrelation time by about three orders of 
magnitude. A more detailed study of the autocorrelations is in preparation. 

To extract the critical exponents we have performed two sets of simula- 
tions; one in the critical region where the correlation length of the infinite 
system is much larger than the linear system size, ^ 3> L, and the other in 
the high-temperature phase, making sure that ^ < L to reduce finite-size 
corrections as much as possible. The quantities we have measured are the 
internal energy, specific heat, correlation length, average cluster size, suscep- 
tibility, and higher order moments. In the critical region this can be done 
most efficiently by combining the single-cluster update with multi-histogram 
sampling techniques, and in the high-temperature phase by using improved 
estimators for measurements. 



2.2 Improved estimators 



An improved "cluster estimator" for the spin-spin correlation function in the 
high-temperature phase, G(x — x') = (s(x) ■ s(x')), is given by |L4 



G(x - x') = 3r—r • s(x) r ■ s(x')Q c (x)Q c (x') 



(5) 



where r is the normal of the mirror plane used in the construction of the 
cluster of size \C\ and Qc{ x ) is its characteristic function (=1 if x G C and 



1 This can be read off from Fig. 2 in ref. 24 1. 
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otherwise). For the Fourier transform, G(k) = J2x G{x) exp( 
implies the improved estimator 



6{k) = W\ 




s(x) cos kx + 



r ■ s(x) sin kx 
\xec J 



-ikx), this 



(6) 



which, for k = 0, reduces to an improved estimator for the susceptibility \ 
in the high-temperature phase, 

\ 2 



G(0) = ~x 




six 



(7) 



Note that we follow in the definition of the susceptibility \ i n t ne high tem- 
perature section the standard convention of omitting the /3-factor. From the 
theoretically expected small k behavior of the inverse Fourier transform, 



G(kY 



E2(l 

Li=l 



cos ki) + (!/£)' 



(8) 



where c is a constant and hi = (27r/L)n$, rij = 1, . . . , L, one may extract the 
correlation length £ by measuring G for a few long-wavelength modes and 
performing least-square fits to (|8]). In our simulations we have measured G(k) 
for n = (0,0,0), (1,0,0), (1,1,0), (1,1,1), (2,0,0), and (2,1,0) (see Sec. 4). It 
is well known that by means of the estimators (0)-(0) a significant reduction 
of variance can only be expected outside the FSS region where the average 
cluster size is small compared with the volume of the system. 



2.3 Histogram techniques 

Even though histogram reweighting techniques have been known since long 
^H , they have gained increasing popularity as a practical tool only quite 



recently . The best performance is achieved near criticality, and in this 



sense the histogram reweighting technique is complementary to the use of 
improved estimators. It is a quite general technique of data analysis based 
on the simple idea of recording whole distribution functions, and not only 
their first few moments (e.g., the average energy and specific heat), as is 
usually done. The energy distribution Vp (E) (normalized to unit area) at 
inverse temperature (3q can be written as 

V Po (E!)=p(E!)e-* E /z(p ), (9) 
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where p(E) is the density of states with energy E. It is then easy to see that 
an expectation value (f(E)) can in principle be calculated for any (3 from 

{mm ~ J? dEPtoWe-W-W ■ (10) 

To keep the notation short we have suppressed the lattice-size dependence of 
Vp (E). 

In practice the continuous energy has to be discretized into bins of size 
AE, and one measures the associated histogram Pp (E) = V/3 (E)AE. Since 
the wings of Pp (E) have large statistical errors, one expects eq. (|H]) (or its 
obvious discrete modification) to give reliable results only for (3 near (3q. If (3q 
is near criticality, the distribution is relatively broad and the method works 
best. In this case reliable estimates from (|T0|) can be expected for (3 values 
in an interval around (3q of width oc L~^ u , i.e., just in the FSS region. The 
Fig. 1(a) shows three typical energy histograms measured at slightly different 
temperatures near the critical point for the lattice size L=48. 

The information stored in Pp Q (E) is not yet sufficient to calculate also 
the magnetic moments (m fc )(/3) with m = \m\ as function of (3 from a single 
simulation at j3 . Conceptually, the simplest way to do so is to record the 
two-dimensional histogram Pp (E, M), where M = mV is the total magneti- 
zation (in general we will denote by small letters quantities per site, and the 
associated total quantities by the corresponding capital letters). Since disk 
space limitations prevented us from doing thatQ we have measured instead 
the "microcanonical averages" 

((m k ))(E) = Y,Pp (E,M)m k /Pp (E), (11) 

M 

where we have used the trivial relation Y^m Pf3o{E, M) = Pp (E). In practice 
this can be done simply by accumulating the measurements of m k in different 
slots or bins according to the energy of the configuration and normalizing at 
the end by the total number of hits of each energy bin. Clearly, once {{m k )) (E) 
is determined, this is a special case of f(E) in eq. (|T0D , so that 

, k) {Q) _ E E ((m k )KE)P, (E)e-^ 



2 As well as from storing the full time series of E and Al, which conceptually would be 
even simpler. 
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Similar to p(E) in theoretically the microcanonical averages ((m k ))(E) 
do not depend on the temperature at which the simulation is performed. 
Due to the limited statistics in the wings of Pp (E), however, there is only 
a finite range around E = (E)((3 ) where one can expect reasonable results 
for ({m k ))(E). Outside of this range it simply can happen (and does happen) 
that there are no events to be averaged. This is illustrated in Fig. 1(b), where 
((m)) (E) is plotted versus energy as obtained from the three runs which gave 
the energy histograms in Fig. 1(a). We see that the function looks smooth 
only in the range where the corresponding energy histogram in Fig. 1(a) has 
enough statistics. 

To take full advantage of the histogram reweighting technique we have 
performed for each L typically three?] simulations at slightly different inverse 



temperatures Using histogram reweighting and jackknife-blocking |]33 
we computed the /3-dependence of the expectation values for all interesting 
thermodynamic observables Oj = 0^(13) plus the associated error AOj. 
To obtain a single expression O = Ol(/3) we then combined the values Oi 
numerically according to the formula 

0=(-^ + -^ + -^)(AO) 2 , (13) 
\(A0 1 ) 2 (A0 2 ) 2 (AO z ) 2 J 



where AO is given by 



1111 



(Aoy (aoo (Ao 2 y (ao 3) 

This expression minimizes the relative error AO/O. The reweighting range 
of each simulation, i.e., that range in which the energy histogram has enough 
statistics to allow for ( |10D to be valid, was determined by the energy values 
at which the histogram had decreased to a third of its maximum value @ . 
From the energy range one can then deduce a corresponding /5-range. Only 
the values inside this /3-window were used for the optimized combination 
(|T3"D, the contributions of the other values were suppressed by giving them 
zero weights. 

We also implemented the optimized histogram combination discussed in 
ref. [^7] to cross-check the results obtained by (|13|) for the specific heat. 

3 Due to technical reasons we had four simulations for the lattices L=24 and L=48, 
compare also Table 1. 
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The Fig. 1(c) shows the weights for the optimal combination of the primary 
energy histograms according to the prescription of ref. ||27|| , i.e., which gives 
the optimal combination of the three estimates for N(E) = p(E)AE. Both 
methods gave comparable results within the statistical errors. We preferred 
our optimization procedure over the optimized histogram addition because 
it is simpler to apply to quantities involving constant energy averages such 
as ((m))(E), and, more importantly, minimizes the error on each observable 
of interest separately. 



3 Results at criticality and finite-size scaling 
analysis 

We investigated in our MC study simple cubic lattices of volume V = L 3 , 
where L=12, 16, 20, 24, 32, 40, and 48. For each L we have made at least 
three simulations at three different temperatures compiled in Table 1. For all 
L we took p = 0.6929, because this is the critical inverse temperature found 



in the recent study of Peczak et al. [23]. We analyzed this run to locate a 



first estimate for the temperatures of the maxima of the specific heat and 
the susceptibility for each L, and used those two temperatures for our other 
two simulations. This choice has the advantage that both locations of the 
maxima scale approximately like L^ 1 ^ which is also the region after which 
the energy distribution Pp{E) has decreased by roughly the same factor for 
all L. We can therefore expect to have always enough overlap to ensure a safe 
reweighting of our three histograms into the /3-region of interest. A further 
advantage is that the two maxima approach T c from different sides. 

For each configuration we recorded the energy histogram Pp(E), using 
90000 bins to discretize the continuous energy range < E < 3V. We 
have checked that this binning is fine enough to ensure negligibly small dis- 
cretization errors. In addition we recorded the microcanonical averages of 
({m))(E), ((m 2 ))(E), and ((m 4 ))(E). These histo grams provided us with all 
necessary information to calculate all thermodynamic quantities of interest. 
Every 10000 measurement steps we recorded a copy of each histogram, so 
that we were able to compute errors by standard jackknife-blocking [[33[ . 
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3.1 Binder parameter Ul{0) 

We have used the histogram reweighting technique to find the /3-dependence 



of the Binder parameter [35 



1-1 



(m 4 ) 
3 (m 2 ) 2 



(15) 



To obtain a single curve Ui{f3) for each lattice size L from the three simula- 
tions at temperatures compiled in Table 1 we used our optimized combination 
of eq. (0). In (|T5| ) we adapt the usual normalization convention, although 
the factor 3 is only really motivated for the Ising (0(1)) model. For general 
0(n) models it is easy to show that in the high-temperature limit Gaussian 
fluctuations around rh = lead to Ul — ► 2(n — l)/3n. For n = 1 this gives 
a zero reference point, but for the Heisenberg model with n = 3 we get the 
quite arbitrary looking limit Ul — > 4/9. In three dimensions, we are for 
low temperatures always in the magnetized phase and hence have for all n 
trivially U L -»• 2/3. 

It is well known [[JIJ that the Ul{@) curves for different L cross around 
(Pc U*) with slopes oc L 1 ^, apart from confluent corrections explaining small 
systematic deviations. This allows an almost unbiased estimate of /3 C , the 
critical exponent u, and U* 
predict f40 | 

U* = 0.59684 . 



Field theoretical expansions in y/e = \/A — D 

(16) 



This follows from the expansion 



(m 4 ) 
(m 2 ) 2 



1 — Xn\/6 



+ Qx z 

12 
B 2 



3T 2 (3/4) 
4T 2 (5/4) 

r(9/4)r(5/4) 



r(7/4)r(3/4) 

1 — x \/6 ^ " 



T(9/4) T(5/4) 

r(7/4) r(3/4) 
Q r 2 (7/4) 
r 2 (5/4) 



g r(7/4) \ 
r(5/4)J 



-4 



(17) 



3 R 



Qx n 



48" 



-R 2 



27 
R 2 



-4 



where x = -1.7650848012^^ = -0.76815456 ^/e, and in the last line 
we have abbreviated R = r(l/4)/r(3/4) = 2.95867512 . . .. Notice also that 
(|15"D can be rewritten as 



2 
3 



a: 



3x 2 



(18) 
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where a 2 = V 2 j3 2 ((m 4 ) — (m 2 ) 2 ) is the variance of the susceptibility. In the 
T — > oo limit we get a 2 , — > |x 2 and at criticality, inserting the field theory 
prediction fllTf), this implies a 2 ~ 0.21x 2 . 

Figure 2 shows the crossings of the Ul{P) curves on a large scale, which 
gives a first estimate of (/3 C , [/*). To extract more precise values of U* and /3 C 
from our data we used that the locations of the crossing point /3 X = 1/T X of 
two different curves Ul(P) and Ul'{P) depend on the scale factor b = L'/L, 
due to the residual corrections to FSS |3^] . To have enough data points for 



a straight line fit we used only the crossing points of the L=12 and L=16 
curves with all the other ones with higher //-value, which gave us 6 and 5 
data points, respectively. The two least-square extrapolations in the plot of 
T x vs 1/ log b shown in Fig. 3 are consistent with each other and gave us the 
values p c = 0.69297(9) and (3 C = 0.69298(13) respectively. Combining the 
two values we obtained 

p c = 0.6930 ± 0.0001. (19) 

The errors on T x were obtained by using the crossings of Ul(P) + AUl((3) 
with Ul'(P) — AUl'(P), where AUl are the errors on Ul obtained via the 
jackknife procedure. Our result ( |19"D is in good agreement with the value 
given in ref. P3] , ft c = 0.6929(1), but significantly higher than estimates 
from analyses of high-temperature series expansions |36 , [37] . 



A similar analysis |35j was used to determine U*, which resulted in U* = 
0.62175(35) and U* = 0.62141(75), repectively, from which we extract the 
final estimate 

U* = 0.6217 ±0.0008. (20) 

The theoretical prediction (fl6| ) based on the -y/e-expansion is about 4% 
smaller than this value. The deviation is somewhat less than for the XY 
model where it is about 6% [[H], 0. 

To determine the derivative dlli/dp, we first used a finite difference ap- 
proximation at our estimate of (3 C . But we observed that the results depended 
very sensitively on the interval of the linear approximation, and on the kind of 
finite difference derivative used, i.e., backward, forward, or symmetric deriva- 



tive. We have therefore chosen another method [|38| to calculate the slope of 



Ul(P) that is less sensitive to systematic errors. We took the thermodynamic 
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derivative of Ul with respect to /3, which can be written as 



W ~ ' { ) ~ W) +{mE) j 

The expectation values for (m k E)((3) at any temperature f3 can be calculated 
from 

Y,E((m k )){E)Pp Q {E)e~^-^ E 

{mtEm = ~ • (22) 

E 

generalizing eq. (|H]). For each L and (3 we obtained in this way three esti- 
mates of ^jg that we combined optimally according to the errors of 

In Fig. 4 we plot ^^(/3 C ), taken at our estimate of (3 C = 0.6930, versus L 
on a log-log scale. From the inverse slope we read off 

v = 0.704 ± 0.006, (23) 

with a quality factor Q=0.61, relying on the linear least-square fit routine 
FIT of ref. ||39|| . For comparison, the field theoretical estimates are v = 
0.705(3) (resummed perturbation series 0]), v = 0.710(7) (resummed e- 
expansion j|2] ) . We noted that the value of v derived in this way was strongly 
dependent on the temperature at which it was extracted. It turned out that 
the systematic errors one gets by choosing a slightly different value for f3 c 
are larger than the statistical errors. The strong /3-dependence of has 
its origin in the large system sizes we used that make the thermodynamic 
quantities vary very rapidly with (3. For comparison, for (3 = 0.6929 (/3 = 
0.6931) we obtained v = 0.696(6) with Q=0.76 (u = 0.712(6) with Q=0.49); 
see also Fig. 8 below. 

3.2 Magnetization and susceptibility 

To extract the ratio of critical exponents (3/v we used that the magnetization 
at (3 C should scale for sufficiently large L like 

(m) oc L- p/v . (24) 
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The slope of the straight line in the log-log plot shown in Fig. 5 of (m) vs L 
at (3 C = 0.6930 gives a value of 

P/v = 0.514 ±0.001, (25) 

with a Q-value of 0.68 of our linear least-square fit routine. We noted again 
that the systematic error of choosing an incorrect (3 C was larger than the 
statistical error coming from the fit. For example, for fits of (m) at the two 
inverse temperatures (3 = 0.6929 and (3 = 0.6931 we obtained (3/v — 0.519(1) 
and (3/v — 0.509(1) with values Q = 0.30 and Q = 0.31, respectively. We 
take the significant lower Q-values as supporting our estimate of (3 C (compare 
also Fig. 8 and the discussion of the Q- values later in this section). Using our 
estimate of v = 0.704(6) we get for the critical exponent of the magnetization 
(3 = 0.362(4). 

On finite lattices an estimator for the magnetic susceptibility per spin is 
given by 

X c = VP((m 2 ) - (m) 2 ), (26) 

where the superscript c stands short for "connected" . In the high-temperature 
phase the true magnetization vanishes, (rh) = 0, and one may also use 

x = V(3(m 2 } for (3 < f3 c . (27) 

Since both expressions should scale at criticality like L 1 ^ = L 2 ~ v a straight 
line fit in a double logarithmic plot of x/L 2 versus L gives an estimate of rj. 
We preferred to visualize r\ as opposed to j/is because r\ is the more sensitive 
parameter. The results for different choices of j3 c are shown in Fig. 6. We 
see that our estimate of f3 c in ( |T9"D seems to be very reasonable, which is 
also supported by the fact that it results in the fit with the highest Q- value. 
The fact that the curves with j3 > (3 C {(3 < (3 C ) bend upwards (downwards) 
is easily understood by recalling that in the low-temperature phase \ scales 
with the volume V = L 3 , while in the high-temperature phase \ eventually 
saturates at a finite value for any L. Close to criticality, as long as the scaling 
variable x = (1 — (3/ (3c)L l / u satisfies |x| 1 and L 1, we expect from FSS 
that the deviations from —77 log L + Co should be simply given by c\x, where 
Co and C\ are approximately constant (apart from small confluent corrections 
oc L~ w ). For all data points in Fig. 6 we have |x| < 0.18 and the linear 
scaling with (1 — /3//3 c ) is obviously satisfied. Moreover, a closer look shows 
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that also the scaling with L 1 ^ is very well confirmed (using our estimate of 
v = 0.704 the deviations from the straight line fit of the curves at (3 7^ (3 C 
for L = 32, 40 and 48 should be stronger than those for L = 24 by factors 
of about 1.5, 2.1 and 2.7, respectively). In fact, in a scaling plot all curves 
would fall on top of each other. The theoretical predictions for 77 are 0.040(3) 
(resummed e-expansion) and 0.033(4) (resummed perturbation series), and 
inserting (p5|) in the scaling relation 77 = 2fi/u — 1 gives 77 = 0.028(2). Our 
direct estimates of 77 for f3 = 0.6930±0.0001 are collected in Table 2. It is also 
interesting to note that 77 gets closer to the theoretical values for (3 < (3 C for 
X, and for {3 > j3 c for x c , hence just in the regions in which both expressions 
are naturally used. Another check on 77 can be obtained from the maxima 
of x°j that should obey the same scaling law as x° and x at (3 C . The results 
are compiled in Table 2. Our three estimates for 77 barely agree in the 2a 
range, but this is mainly due to the low estimate of 77 coming from the x c ~ 
fit. This behavior of the x° scaling law is in agreement with the findings of 
refs. p3| , |32| , who also observed that the x°-fit results in a noticeably lower 



77 value. As our final estimate we take the result from the best fit of x at 
P c = 0.6930, yielding 

77 = 0.0271 ± 0.0017. (28) 

This in turn implies 7/V = 2 — 7/ = 1.9729(17) and, using again our value of 
v, 7 = 1.389(14). 

From the temperature locations of the (connected) susceptibility 

maxima one can get another estimate of (3 C by using the scaling prediction 

T x o max = T c + aL~ l l u + Using our previously determined value of v = 0.704 

we get from the linear fit shown in Fig. 7 the estimate (3 C = 0.6930(3), with 
Q=l, which is in excellent agreement with our result ( |i~9|) for (3 C obtained 
from the Ul crossing points. The data for the fit is compiled in Table 3. 

When we did our least-square fitting of quantities extracted at a par- 
ticular (3, we noticed that the Q-value can also be used as a very sensitive 
measure to obtain an estimate of the critical temperature, whose error is 
comparable to the one obtained from the fits to the Ul crossing points. Ob- 
viously if one extracts data at temperatures away from (3 C the subleading 
corrections to the scaling laws become more important and worsen the fit 
substantially, especially if large lattices are used. In Fig. 8 we plot the qual- 
ity factor Q obtained from log-log linear least-square fits of the quantities 
^jp-{f3)-, {m)((3),x{f3), and x c (P) vs the lattice size L, extracted at different 
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inverse temperatures. The most striking results come from the Q(/5)-curves 
of (m) and Xi that seem to have the most severe f3- dependence. Their peak 
locations are exactly at (3 — 0.6930 with an estimated error of ±0.0001, 
whereas the Q(/3)-curves of ^j^(P) and x c (P) do not peak so sharply at some 
(3 and are asymmetric with respect to their peak location, falling off very 
slowly on the high-temperature side. The slow variation of the Q-value of 
the ^#(/3)-fits indicates that Binder's method to determine v is less sensitive 
to the choice of (3 C than the FSS method for (3/v and j/u, as is expected on 



theoretical grounds ||35| . This is reflected in our data for u, where the sys- 
tematic errors for different choices of the critical coupling are of the order of 
the statistical errors, whereas for they are five times as big. 



3.3 Energy and specific heat 

Assuming hyperscaling, a = 2 — Du, to be valid and inserting our estimate 
v = 0.704(6), we get a negative value for the critical exponent of the spe- 
cific heat, a = —0.112(18). This implies a cusp like singularity at T c , but 
no divergence at all. Although at first sight this sounds numerically very 
convenient, it turns out that the specific heat is the hardest observable to 
analyze. On one hand it is quite easy to get the peak height with good preci- 
sion due to the smoothness of the specific heat. On the other hand, however, 
the scaling behavior is quite unclear for a < 0, since one has to deal with 
large, non-universal background terms. For the peak location, Tc max (L), the 
scaling behavior is of the usual form T<7 max (L) « Tc + aL~ x l v ', but now the 
smoothness of the peaks makes it numerically difficult to determine Tc max (L) 
with high precision. 

We calculated the specific heat in two different ways, once from the energy 
fluctuations as 

C = V73 2 «e 2 >-(e> 2 ), (29) 

and the other time as a finite difference derivative approximating C = de/dT. 
Both definitions resulted in similar curves, and this time the choice of the 
finite difference scheme did not affect the results. Because the curves of 
definition (|29"D were smoother in appearance, we decided to use solely the 
first definition to extract our scaling information. The maximum of the 
specific heat should scale as 

C max (L) « CSL - aL«/», (30) 
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where C^| x is a regular background term. Using the routine MRQMIN of 
ref. [3{| we obtained from a three-parameter fit the values C^| x = 4.17(97), 
a = 3.98(17), a/v = -0.33(22), with a Q-value of Q=0.69. Using our 
estimate of v — 0.704(6), this results in a = —0.23(16). This estimate is 
clearly larger than the one obtained through hyperscaling but, due to its 
large error bar, it is still consistent. By imposing the hyperscaling value of 
a/v = 2/v — 3 = —0.159, we also did a linear fit of C max (L) vs L a / U that 
gave us C^ e f x = 5.79(12), and a = 5.00(19) withfl Q = 0.71. Both curves fit 
the data almost equally good as can be inspected in Fig. 9. 

The fit of 7c max (L) versus L' 1 ^ is shown is Fig. 7 . Assuming v = 0.704, 
our linear fit routine gave (3 C = 0.6925(9) with Q=0.80. Due to the large 
error bars on Tc max this estimate has the largest statistical scatter, but is 
still in agreement with our previous estimates of (3 C . The data for the fit can 
be found in Table 3. 



4 Results in the high-temperature phase 

Let us now turn to the analysis of our second set of simulations which were 
performed in the high-temperature phase. Here we have mainly concentrated 
on measurements of the spatial correlation length, £, and the susceptibility, 
X, using the improved estimators of Sec. 2.2. From least-square fits to the 
critical divergencies of £ and \ as function of T or (3 we can get independent 
estimates of the critical coupling (3 C and of the critical exponents v and 7. 
This of course requires to keep finite-size corrections as small as possible. We 
therefore have chosen the linear size of the lattices to satisfy the condition 
L > 8£, which proved in related studies [[0], to be a safe condition. In a 
few cases we have checked this once again by performing test runs at smaller 
L ~ 5 — 7£. As a result, to avoid the most severe finite-size corrections, also 
for this model it would probably be sufficient to choose L pa 6£. 

Our final data set consists of 18 points between (3 = 0.650 and (3 = 0.686 
(corresponding to correlation lengths £ = 3, . . . , 12, see below), using lattices 
of sizes between 32 3 and 100 3 . The simulation parameters together with the 
statistics and the results are compiled in Table 4. The statistics is given in 
terms of N, the number of measurements of the standard estimators which, 

4 Although the Q-value for the 2-parameter fit is better, its total % 2 (=2.9) is larger 
than that of the 3-parameter fit (=2.2). 
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on the average, were taken after / x V spins are flipped. The number of 
Metropolis equivalent sweeps is then simply / x N. The average cluster size 
(|C|) is given in terms of the ratio (]C|)/xi mp which is only very slowly varying 
over a wide temperature range. The column labeled with C shows the specific 
heat calculated from the energy fluctuations. For comparison we give for the 
susceptibility the averages over standard as well as over improved estimators. 
We see that the errors on x~ imp are smaller by about a factor of 2 — 3. It should 
be mentioned that the errors on x could be somewhat reduced by doing 
measurements more frequently (i.e., by choosing smaller values of /). The 
reason is that the (integrated) autocorrelation time for \ is still very small 
on the time scale at which the measurements are taken 0.8, translating in 
Metropolis equivalent units to 0.8 x / ~ 0.1). On the other hand it should 
be kept in mind that each standard measurement takes 0(V) operations, so 
that too small a factor / would slow down the simulation considerably. The 
correlation length in the last column was extracted from fits to the inverse 
Fourier transform of the spatial correlation function as described in Sec. 2.2. 
We always used the lattice momentum squared, X^ =1 2(l — cos/cj), ki = 
(2n/L)ni, as independent variable. In Fig. 10 this procedure is illustrated 
for our largest lattice of size 100 3 . We have checked that within error bars 
the estimates for £j mp do not depend on how many Fourier coefficients we use 
in the fits. This observation is also supported by the very reasonable values 
of the goodness-of-fit parameter for all fits. Notice that even the simplest 
expression, £ imp = (G(0)/G(1) — l) 1 / 2 /2 sin(7r/L) (involving no fit at all), 
can be used. Our final results for £ imp in Table 4 are from fits using all six 
possible k values up to k = (2ir/L)(2, 1,0), since the error estimates trivially 
decrease with the number of points taken into account. 

In what follows we describe the least-square fit techniques employed to 
determine the critical coupling and the exponents v and 7 from the raw data 
in Table 4. Although we are quite confident in the accuracy of the correlation 
length results, we shall first consider the susceptibility data which, a priori, 
is more reliable since no intermediate analyses are involved. 

4.1 Susceptibility 

Starting with the simplest ansatz, 

x(P) = Xoa-P/Pc)~\ (31) 
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we obtained from a non-linear three-parameter least-square fit to all 18 data 
points (using x imp ) 



p c = 0.69294 ±0.00003, (32) 
7 = 1.391 ±0.003, (33) 
Xo = 0.955 ±0.006, (34) 

with a chi-squared of x 2 — 7.83, corresponding to a goodness-of-fit param- 



eter Q = 0.93, which clearly justifies a posteriori the ansatz (|31|) . The 
very precise estimate for (3 C is only little lower than our previous value (Tl9|), 
P c = 0.6930(1), derived from FSS analyses of the Ul crossing points. The 
exponent 7 agrees very well with the field theoretical estimates, 7 = 1.386(4) 



(resummed perturbation series [41]), 7 = 1.390(10) (resummed e-expansion 
[fm). The amplitude Xo is close to the mean-field value, x{f F = 1, and agrees 
surprisingly well with old analyses of 9-term high-temperature series (HTS) 
expansions by Ritchie and Fisher f36j (RF) which gave x^ F = 0.96647(5) 
(using (3f F = 0.6916(2) and 7 RF = 1.38(2) as input). We could not trace 
any more recent published number to compare with. From the HTS analy- 



ses in ref. [37 , however, it is straightforward to derive another comparative 
value for xo- Expanding the ansatz for the dominant singularity (j3~lD into 
a power series in /3, inserting the estimates for the critical parameters given 
in Table 4 of ref. |7|, p c = 0.6924(2), 7 = 1.387(4) (Pade approximants) 
and p c = 0.6925(1), 7 = 1.395(5) (ratio method), and comparing with the 
coefficients of the HTS expansion |37], we obtain a sequence of amplitude 
values that stabilize with increasing order of the expansion at Xo TS ~ 0.9508 
and Xq TS ~ 0.9326, respectively. Moreover, using our values ([32]) and ( |33|) 
for P c and 7, we obtain by the same procedure a very smooth sequence for 
Xo yielding asymptotically xo — 0.9501. The smoothness indicates that the 



simplest ansatz ( 31|) together with our estimates (|32|)-(|341) represents an ex- 



cellent extrapolation of the HTS expansion, whose quality is comparable (if 



not better) with the results given in ref. [[37] . 

As a further self-consistency check of the range of validity of the ansatz 
( |31~]) we tested for the asymptotic scaling region by successively discarding 
more and more data points with small correlation length. As is demonstrated 
in Fig. 11 for P c , only a weak downward trend is observable which, compared 
to the error bars, is hardly significant. We can thus conclude that the ansatz 
( pT]) passes all usual self-consistency checks and that the estimates (p2|)-(p^) 
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should thus be reliable. 

It is of course surprising that the simple ansatz (|3T|) works even for data 
points with correlation lengths as small as £ ~ 3. Theoretically we would have 
expected to see confluent correction terms of the type Xconf(l — /^/A:)~ 7+Al , 
where A x = ujv ~ 0.55 0Qis the confluent correction exponent [Q, and 
analytic correction terms of the Darboux type, Xanai(l — PI Pc)~~ <+1 , leading 
to the more general ansatz 

X(P) = Xo(l - P/Pc)^ + Xconf (1 " /3/A)" 7+Al + Xanal(l " P / Pc)^ ■ (35) 

Including both correction terms as free parameters into the fit routine is a 
hopeless enterprise. Instead, we first tried fits with p c held fixed at values 
around 0.69290 in steps of 0.00001, and 7, Xo, Xconf, and Xanai as free parame- 
ters. The quality of the fits remained high for a large range of /5-values. The 
best fit with Q = 0.97 was obtained for (3 C = 0.69281 and gave 7 = 1.365(25), 
Xo = 1.06(17), Xconf = -0.11(79), and Xanai = -0.25(1.18). We see that the 
correction terms come out to be consistent with zero, but that the error bars 
are fairly large due to the increased number of free parameters, so that no 
definite conclusion can be drawn from these fits. 

We therefore decided to perform two further types of fits. First, we added 
to the leading term in (^) only the analytic correction term (i.e., we enforced 
the amplitude of the confluent correction term to be zero), and second, we 
tried to use only the confluent correction term with fixed Ai = 0.55 (i.e., 
we put the amplitude of the analytic correction equal to zero). Fitting thus 
all 18 data points to the ansatz (|35| ) with either Xconf or Xanai held fixed at 
zero we get for the other amplitude Xanai = —0.40(27) or Xconf = —0.35(26), 
respectively. The Q-value improves only marginally, and the amplitudes are 
still almost consistent with zero. The other parameters change slightly, but 
due to the much larger error bars than for the simple fit (^1|) they are still 
compatible. 

The latter fits are non-linear four-parameter fits which are usually quite 
difficult to stabilize. To cope with this problem we followed ref. [TTJ and 



used the following method. For any pair of values P c , 7 we first minimized in 
the linear parameters Xo and Xanai (or Xconf), which can be done exactly (up 
to round-off errors). This yields a chi-squared function x 2 (Pc, 7, Xo(Pc 7), 
Xanai(/? c , 7)) = X 2 (Pc,l) which depends on two non-linear parameters only 

5 The precise values are u = 0.78(2), A 1 = 0.55(3) 0; u) = 0.79(4), Ai = 0.56(4) (42). 
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and can be minimized reliably by standard subroutines. Finally we used 
the so determined estimates of (3 C , 7, Xo, Xanai (or Xconf) as initialization of 
the non-linear fit-routine MRQMIN of ref. |3"9]] , which yields then slightly 
improved parameter values and standard error estimates. The fits with 
7,XcbXconf an d Xanai as free parameters were performed similarly. 
One can also rewrite eq. (|35|) as a function of T, 



X(T) = Xo(T/T c - ir + xLf (T/T e - ir +Al + ^(r/r e - l)"^ 1 , (36) 

with Xo = Xo, Xconf = Xconf, and x' ana i = Xanai + 7Xo- This form is obtained 
from (|35| ) by a simple change of variable f3 — > 1/T, expanding around T c , 
and keeping the same correction terms. We believe that (|36| ) with Xconf = 
Xanai = is, a priori, as justified as the simplest /3-dependent ansatz (]31f) . 
Since HTS analyses are based on expansions in j3, there is definitely a bias 
to prefer the ansatz (|31|). We are not aware, however, of a mathematically 
sound justification for this choice. For what we think is a counter-example 
see our correlation length analysis below. 



Because the simplest /3-dependent ansatz ( fHf ) works so good, it is thus 
clear that for the temperature dependent fit ansatz we have to take into 
account the analytic correction term. As anticipated, trying to fit the sim- 
plest ansatz (^) with Xconf = Xanai = to all 18 data points results in a 
fit with a comparatively poor chi-squared of x 2 = 21.66, corresponding to 
Q = 0.12. By again successively discarding the data points with small cor- 
relation length, the quality of the fit improves rapidly, and in Fig. 11 we see 
a weak trend that the values of j3 c from the two types of fits may approach 
each other asymptotically. 

The fit to all 18 data points with x'anai as a f ree parameter agrees perfectly 
with the corresponding /5-dependent ansatz (also with Xanai as a free param- 
eter), yielding the expected correction amplitude x ana i = 1-03(26). The fit 
with Xconf as f ree parameter is only slightly worse. The value of Xconf = 0-90 
(instead of the expected result Xconf = X«mf ~ 0) indicates, however, that 
in this case the confluent correction term trys to account for the suppressed 
analytic correction term, and that it can indeed mimic the analytic behavior 
very well. 

The main result of all our efforts is that on the basis of standard statistical 
tests, the simple ansatz (|3lD turns out to be justified. We emphasize this 



point because a priori it is not at all clear that the confluent correction is 
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extremely small and that also the analytic correction is negligible if the proper 
variable {(3 in the case of x) is chosen. Although we feel that the corrections 
cannot completely be ignored, one would need more accurate data to account 
for these corrections in a reliable way. 



4.2 Correlation length 

We have repeated the same type of analysis for the correlation length data. 
Starting again with the simplest possible ansatz taking into account only the 
leading singularity, we find this time that the temperature dependent ansatz, 

£(T) = &(T/T C - I)"", (37) 

is well behaved. It should be emphasized that, as far as the role of (3 and T 
is concerned, this is just opposite to the situation for the susceptibility. The 
fit to all 18 data points in Table 4 yields 

(3 C = 0.69281 ±0.00004, (38) 
v = 0.698 ±0.002, (39) 
Co = 0.484 ±0.002 (40) 

with a chi-squared of \ 2 = 8.14, corresponding to a goodness-of-fit parameter 
Q = 0.92, which may again be taken as a posteriori justification of the simple 
ansatz fl3~7|). The estimate for (3 C is somewhat smaller than the previous 
estimates, and also the exponent v is only barely compatible with our FSS 
value v = 0.704(6) and with the field theoretical estimates, v = 0.705(3) 



(resummed perturbation series [41]), v = 0.710(7) (resummed e-expansion 
[0), but it still lies within the 2a error interval of these estimates. If we 
extract £ from fits through the four lowest k values only, we obtain the same 
central values as in (|38|)-(|40|) with about 1.5 times bigger error bars and a 
Q- value of 0.99. In Fig. 11 we see that fits with the ansatz (|37|) are very 
stable against discarding more and more data points with small correlation 
length. 

Correspondingly, if we include into the T-dependent ansatz one cor- 
rection term at a time, as described in the susceptibility analysis, we ob- 
tain amplitudes that are fully consistent with zero, ££ onf = —0.048(65) and 
£anai = —0.065(76). For the corresponding /^-dependent fit we thus ex- 
pect £ con f ~ and £ ana i = <^ nal — z/£q = —0.267. The latter value is ex- 
tremely well reproduced by the /5-dependent fit with analytic correction term 
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(£anai = —0.272(76)). Similar to the susceptibility the confluent amplitude, 
however, comes out much too large in an attempt to mimic the analytic 
correction term. We can thus conclude that also for the correlation length 
confluent corrections are negligible. Analytic correction terms are important 
in the /3-dependent ansatz but negligible in the T-dependent ansatz, just 
opposite to the situation for the susceptibility. 

It is noteworthy that, using correction terms, the value of v slightly in- 
creases, just opposite to the case of the susceptibility data, where the value 
of 7 decreases when corrections are included. It is impossible to find an ob- 
jective criterion for the systematic errors, but, with what we experienced, we 
feel that they are at least of the order of the statistical errors, in particular 
for the susceptibility, where the correction terms seem to be somewhat more 
important than for the correlation length. 

Combining the estimates for (3 C from the susceptibility and correlation 
length data in the high-temperature phase we obtain a final estimate of 

(3 C = 0.69288 ± 0.00004, (41) 

which is slightly smaller than our crossing value of 0.6930(1), but in agree- 



ment with the MC estimate by Peczak et al. [23 . 



4.3 The exponent 77 

Finally, we have combined the scaling behavior of £ and \ to get a direct 
estimate of the exponent 77 = 2 — 7/1/ from the relation \ °^ £ '> or 

log(x/£ 2 ) =c-77log£, (42) 

where c is a constant. By plotting x/£ 2 vs £ on a log-log scale as in Fig. 12(a) 
we thus expect asymptotically for large £ a straight line with slope — r\. We see 
that up to £ ~ 12 our data still give a curved line, indicating that corrections 
to asymptotic scaling cannot yet be neglected. Taking a linear envelope to 
the last few points as rough estimate, we obtain rj « 0.05. Alternatively we 
can define an effective exponent = - log (xi+i^/Mf+i) I log as 
the local slope between two points in Fig. 12(a). Since this gives quite noisy 
results we have actually computed the local slopes from linear fits through 
three and four points, denoted by 77® and 77^ , respectively. These effective 
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exponents are plotted vs £ on a logarithmic scale in Fig. 12(b). Again we see 
that up to correlation lengths of the order of £ « 12 the effective exponents 
have clearly not yet reached the asymptotic value (= if). However, assuming 
a monotonic behavior as £ — > oo we get at least an upper bound, 77 < 0.05. 
Recall the field theory values which are rj = 0.033(4) (resummed perturbation 
series r] = 0.040(3) (resummed e-expansion [|42"fp|). 

Estimates for the other critical exponents can be derived using (hyper-) 
scaling laws. We obtain (5 = \v-\l = 0.352(4), and a = 2-3u = -0.094(6) 
(compare also Table 5). 



5 Concluding remarks 

In conclusion we have shown that the smg/e-cluster update eliminates critical 
slowing down for the three-dimensional Heisenberg model almost completely. 
As for the 3D Ising and XY model we expect that it is more efficient than 
the multiple-cluster algorithm, but this has not yet been explicitly verified 
due to the lack of data for the multiple-cluster algorithm. Combined with 
histogram reweighting and optimization techniques, finite-size scaling analy- 
ses allow a precise Monte Carlo determination of the critical exponents of the 
3D Heisenberg model, whose accuracy is comparable with the best estimates 
coming from field theoretical methods. Direct analyses of thermodynamic 
measurements based on improved estimators in the high-temperature phase 
yield compatible results. For the reader's convenience we have compiled 
some relevant sources on the critical parameters of the classical 3D Heisen- 
berg model in Table 5. Overall, our results are in good agreement with the 
Monte Carlo values reported recently by Peczak et al. [p3| |. In particular we 



confirm that the value for f3 c on a simple cubic lattice is significantly higher 
than previous estimates coming from analyses of high-temperature series ex- 
pansions. We would like to remark, however, that a very recent analysis 
4^fl of extended series expansions |48| (up to 14th order), using more refined 



Pade-approximant techniques, is consistent with the MC estimates. 



6 In a recent recalculation |J45[ of all Feynman graphs contributing to the highest order 
of the e expansion ( C(e 5 ) ) several errors were corrected. A subsequent reanalysis [Q 
of the resummed series gave a slightly smaller value for rj. Compared with the error bar, 
however, this change is negligible. 
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Tables 



L 


P 

1 


N 


f 


( 


C 


) 


T 


12 


0.6783 


1430249 


0.045 




79 




0.8 


12 


0.6929 


5222274 


0.071 


122 


1.2 


12 


0.7009 


931277 


0.086 


149 


1.4 


16 


0.6872 


242804 


0.122 


167 


1.4 


16 


0.6929 


217393 


0.157 


213 


1.8 


16 


0.6953 


206676 


0.174 


241 


1.9 


20 


0.6862 


270067 


0.110 


220 


1.1 


20 


0.6929 


332725 


0.208 


332 


1.8 


20 


0.6947 


284544 


0.184 


369 


1.8 


24 


0.6872 


709711 


0.130 


300 


1.0 


24 


0.6929 


472831 


0.207 


476 


1.6 


24 


0.6953 


350379 


0.244 


563 


1.7 


24 


0.6972 


355572 


0.274 


631 


1.7 


32 


0.6881 


356881 


0.126 


460 


1.1 


32 


0.6929 


351394 


0.206 


842 


1.8 


32 


0.6965 


360726 


0.215 


1174 


1.9 


40 


0.6904 


255563 


0.106 


851 


1.1 


40 


0.6929 


216251 


0.163 


1304 


1.5 


40 


0.6953 


232234 


0.194 


1774 


1.8 


48 


0.6914 


131620 


0.145 


1339 


1.2 


48 


0.6929 


138729 


0.200 


1840 


1.6 


48 


0.6941 


126643 


0.165 


2287 


1.6 


48 


0.6968 


174441 


0.178 


3287 


1.7 



Table 1: Measurement statistics: N is the number of measurements, taken 
after / x L 3 spins are flipped (on the average), (\C\) is the average cluster 
size that is measured after each update step, and r is a rough estimate of 
the autocorrelation time in Metropolis-equivalent units. 
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X 


x c 




Xmax 


p 


V 


Q 


V 


Q 




V 


Q 


0.6929 


0.0364(17) 


0.36 


0.0086(44) 


0.71 




0.0231(61) 


0.30 


0.6930 


0.0271(17) 


0.78 


0.0156(44) 


0.69 








0.6931 


0.0178(17) 


0.43 


0.0237(44) 


0.48 









Table 2: Results for rj obtained from log-log fits of xjL 2 and x c jL 2 versus 
the lattice size L at three temperatures near (3 C . Included is also the estimate 
of r) coming from a fit of the maximum x m ax °f the connected susceptibility 
(see text). Q denotes the quality factor of the least-square fit routine FIT of 
ref. ||. 



L 


T-yC 

A. max 


Xmax 


*— 'max 


Cmax 


12 


1.4747(35) 


5.260(38) 


1.4243(101) 


2.426(11) 


16 


1.4629(28) 


9.242(100) 


1.4362(45) 


2.579(26) 


20 


1.4579(13) 


14.25(10) 


1.4393(39) 


2.711(20) 


24 


1.4550(11) 


20.80(09) 


1.4355(37) 


2.778(15) 


32 


1.4505(11) 


36.53(22) 


1.4386(35) 


2.936(34) 


40 


1.4486(08) 


56.45(39) 


1.4430(21) 


2.959(46) 


48 


1.4474(06) 


81.48(63) 


1.4415(15) 


3.090(50) 



Table 3: Data for the temperature locations and the absolute values of the 
maxima of the magnetic susceptibility and the specific heat. The data was 
extracted from the optimally combined curve of the three runs for each lattice 
size (compare text). 
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L 


T It. 
( Simp 


iV/10 3 


f 


(|C|)/Y imD 
\ 1 1 / / /V. imp 


C 


Y 

A 


Ximp 


•Cimp 


0.650 


32 


9.9 


79.5 


0.270 


0.7556 


1.074(13) 


46.00(15) 


45.697(50) 


3.2345(26) 


0.655 


32 


9.0 


263.6 


0.134 


0.7547 


1.137(11) 


54.43(12) 


54.280(53) 


3.5412(24) 


0.660 


32 


8.1 


93.0 


0.271 


0.7539 


1 222(13) 


65 95(21) 


66 095(89) 


3 9336f34) 


0.665 


40 


9.0 


400.0 


0.133 


0.7531 


1.310(11) 


83.07(14) 


83.155(65) 


4.4343(24) 


0.670 


40 


7.8 


138.8 


0.134 


0.7523 


1.382(19) 


109.32(33) 


109.59(18) 


5.1260(53) 


0.673 


50 


8.8 


89.5 


0.133 


0.7519 


1.453(23) 


132.77(47) 


132.76(22) 


5.6676(62) 


0.675 


50 


8.2 


114.4 


0.132 


0.7517 


1.470(24) 


153.92(48) 


153.99(26) 


6.1197(64) 


0.676 


60 


9.4 


132.3 


0.133 


0.7515 


1.539(21) 


167.07(47) 


166.69(21) 


6.3740(55) 


0.677 


60 


9.0 


68.5 


0.133 


0.7515 


1.560(29) 


181.71(72) 


181.16(34) 


6.6584(83) 


0.678 


60 


8.6 


170.5 


0.133 


0.7513 


1.562(20) 


198.52(49) 


198.74(25) 


6.9869(58) 


0.679 


60 


8.2 


156.3 


0.133 


0.7512 


1.619(21) 


218.82(58) 


218.46(32) 


7.3333(68) 


0.680 


60 


7.8 


95.3 


0.162 


0.7511 


1.700(28) 


242.04(79) 


242.17(44) 


7.7359(87) 


0.681 


70 


8.5 


295.5 


0.133 


0.7510 


1.710(17) 


270.68(51) 


270.96(27) 


8.2013(52) 


0.682 


70 


8.0 


154.9 


0.133 


0.7508 


1.776(23) 


304.38(79) 


305.91(46) 


8.7244(81) 


0.683 


80 


8.6 


73.4 


0.157 


0.7508 


1.765(31) 


348.5(1.3) 


349.26(63) 


9.346(11) 


0.684 


80 


7.9 


138.4 


0.133 


0.7506 


1.778(25) 


405.9(1.2) 


405.92(62) 


10.0988(94) 


0.685 


90 


8.2 


174.0 


0.133 


0.7505 


1.903(24) 


477.9(1.2) 


478.60(64) 


10.9857(92) 


0.686 


100 


8.3 


77.4 


0.133 


0.7504 


1.860(36) 


575.0(2.0) 


576.9(1.1) 


12.093(15) 



Table 4: Results in the high-temperature phase. iV is the number of measure- 
ments of the non-improved observables, taken after / x L 3 spins are flipped 
(on the average). C denotes the specific heat and (|C|) is the average cluster 
size. The cluster size and the improved observables are measured after each 
cluster- up date step. 
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field theory 


MC simulations 




critical 


g-expansion 


e-expansion 


Metropolis MC 


this 


study 


parameter 


ref. [41] 


refs. 40, 42 1 


ref. H 


FSS 


HT 


Pc 






0.6929(1) 


0.6930(1) 


0.69288(4) 


U* 




0.59684 


0.622(1) 


0.6217(8) 




V 


0.705(3) 


0.710(7) 


0.706(9) 


0.704(6) 


0.698(2) 


a 


-0.115(9) 


-0.130(21) 


-0.118(18) 


-0.112(18) 


-0.094(6) 


ajv 


-0.163(12) 


-0.183(28) 


-0.167(36) 


-0.159(24) 


-0.135(9) 




0.3645(25) 


0.368(4) 


0.364(7) 


0.362(4) 


0.352(4) 


P/v 


0.517(6) 


0.518(11) 


0.516(3) 


0.514(1) 


0.504(5) 


7 


1.386(4) 


1.390(10) 


1.390(23) 


1.389(14) 


1.391(3) 


rj 


0.033(4) 


0.040(3) 


0.031(7) 


0.027(2) 


< 0.05 



Table 5: Various sources for estimates of the critical parameters for the 
classical 3D Heisenberg model. For the MC simulations scaling relations 
were used to obtain estimates for a (all) and (3 (only HT), and the FSS 
values of the exponents 7 and (3 were calculated from the measured ratios, 
using the estimate for v. The field theory estimates of ratios with v are 
calculated from the values and errors of the critical exponents. 
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Figure Headings 



Fig. 1: (a) The energy histograms for L = 48 at the three simulation tem- 
peratures, (b) The constant energy averages ((m))(E) as computed 
from the three simulations yielding the histograms in (a), (c) The 
weights for the optimal combination of the histograms according to the 
procedure in ref. |27| . 

Fig. 2: The Binder parameter Ul vs (3. The values of Ul{@) were obtained 
by reweighting and optimized combining the results of our three simula- 
tions at different temperatures for each lattice size L. The simulations 
were performed at f3 = 0.6929 (the critical inverse temperature found 
by Peczak et al. |||), and at the positions of the maxima of the specific 



heat C and the susceptibility x c ; respectively. 

Fig. 3: Estimates for T c , coming from plotting the T-values of the crossing 
points of the L—12 and L=16 curves of the Binder parameter Ul vs 
the inverse logarithm of the scale factor b = L'/L. The extrapolation 
leads to an estimate of (3 C = 0.6930(1). 

Fig. 4: The thermodynamic derivative calculated at /5 C =0.6930 vs lat- 
tice size L in a double logarithmic plot. The slope of the linear least- 
square fit gives an estimate for the critical exponent of the correlation 
length, v = 0.704(6), with a quality Q=0.61. 

Fig. 5: Double logarithmic plot of the magnetization (m) at (3 C = 0.6930 
vs lattice size L. The values of (m) were obtained by reweighting and 
optimized combining (see text) of our three simulations at different 
temperatures for each lattice size L. The slope gives an estimate for 
P/u. From the fit we obtain /3/u = 0.514(1), with quality Q = 0.68. 

Fig. 6: log(x/L 2 ) vs log(L) at inverse temperatures /3=0.6925 (□), /3=0.6927 
(O), ,3=0.6929 (O), ,3=0.6930 (x), ,3=0.6931 (A), /?=0.6933 (A), and 
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,3=0.6935 (o). The straight line corresponds to a linear least-square fit 
at ,3=0.6930, that gives 77=0.0271(17) with Q=0.78. 

Fig. 7: Variation of the pseudo transition temperatures T x ^ (L) and 7c max (L) 
with L^ 1 ^, where z/ = 0.704(6) is our estimate obtained in Fig. 4. The 
fits yield estimates of (5 C = 0.6930(3) (Q = 1.0) and & = 0.6925(9) 
(Q = 0.80), respectively. 

Fig. 8: The Q- value (quality-factor) of linear least-square log-log fits of the 
quantities X C (P), ( m ){l3), and ^f((3) vs the lattice size L, ex- 

tracted at different /5-values in intervals of A(3=0. 00002. 

Fig. 9: Lattice-size dependence of the maxima of the specific heat, C max - 
The solid curve is a non-linear three-parameter fit, whereas the dashed 
curve comes from a linear fit with fixed a/u = 2/u — D = —0.159, 
employing hyperscaling arguments and our estimate of v — 0.704(6). 

Fig. 10: Fit of G(k)- 1 to c[E- =i 2(l-cosfc j ) + (I/O 2 ] on a 100 3 lattice used 
to compute the correlation length £. 

Fig. 11: The critical coupling f3 c as estimated from non-linear three-para- 
meter fits to £ and x assuming the leading power-law singularity only, 
written as function of (3 or T. The x-axis shows how many data point 
are taken into account in these fits (i.e., from right to left more and 
more points with small correlation length are discarded). In order to 
disentangle the error bars two curves are slightly displaced in x. 

Fig. 12: (a) Log-log plot of x/£ 2 vs £. The slope of the linear envelope 
for large £ gives 77 ~ 0.05. (b) The effective critical exponent 77^ vs 
£ on a logarithmic scale, using two different discretization schemes. 
The asymptotic value of i] cS as £ — > 00 is an estimate for the critical 
exponent r\ = 2 — 7/1/ . 
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